% File: softdec.m
% ---------------
% This script implements soft decision process.

% M_mod  = 4; % QPSK
% nsymb  = 2048;

SNR    = snr_in(:);
MSE    = 1./SNR;
M_LIM  = length(MSE);
SER    = zeros(size(MSE));
VAR_SOFT  = zeros(size(snr_in));
VAR_SOFT1 = zeros(size(VAR_SOFT));
% VAR_SOFT2 = zeros(size(VAR_SOFT));
numerr    = zeros(length(SNR)); % Compute error rate for each subframe
h_waitbar = waitbar(0,'Please wait...');
m = 1;
while(m<=M_LIM)
%     numTrials = 10000;
    var_soft  = zeros(numTrials,1);
    var_soft1 = zeros(size(var_soft));
    var_soft2 = zeros(size(var_soft));
    
    n = 1;
    while (n<=numTrials)
        x       = randi([0,M_mod-1], nsymb, 1);
        x_pload = modulate(modem.qammod(M_mod),x);
        xtxPower_mean    = mean(abs(x_pload).^2);
        xtxPower_mean_dB = 10*log10(xtxPower_mean);
        noisePower_dB    = xtxPower_mean_dB-SNR(m);
        noisePower       = 10^(noisePower_dB/10);
        noisePow_norm    = noisePower/xtxPower_mean;
        s_pre_det        = awgn(x_pload, SNR(m), xtxPower_mean_dB);
        
%         s_pre_det_real = real(s_pre_det);
%         s_pre_det_real_soft_pos = exp(-1*(s_pre_det_real-1).^2/(MSE(m)*xtxPower_mean));
%         s_pre_det_real_soft_neg = exp(-1*(s_pre_det_real+1).^2/(MSE(m)*xtxPower_mean));
%         s_mean_real    = (s_pre_det_real_soft_pos - s_pre_det_real_soft_neg)./(s_pre_det_real_soft_pos + s_pre_det_real_soft_neg);
%         s_pre_det_imag = imag(s_pre_det);
%         s_pre_det_imag_soft_pos = exp(-1*(s_pre_det_imag-1).^2/(MSE(m)*xtxPower_mean));
%         s_pre_det_imag_soft_neg = exp(-1*(s_pre_det_imag+1).^2/(MSE(m)*xtxPower_mean));
%         s_mean_imag    = (s_pre_det_imag_soft_pos - s_pre_det_imag_soft_neg)./(s_pre_det_imag_soft_pos + s_pre_det_imag_soft_neg);
%         s_mean_test    = s_mean_real + 1i*s_mean_imag;
%         real_mean      = mean(s_mean_real);
%         real_var       = mean(s_mean_real.^2 - real_mean.^2);
%         imag_mean      = mean(s_mean_imag);
%         imag_var       = mean(s_mean_imag.^2 - imag_mean.^2);
%         s_mean_var     = real_var + imag_var;
%         
%         var_soft(n)  = mean(abs(s_mean_test - x_pload).^2)/xtxPower_mean;
%         var_soft1(n) = mean(2 - abs(s_mean_test).^2)/xtxPower_mean;
%         var_soft2(n) = 1 - s_mean_var/xtxPower_mean;

% -----------------------------------------------------------------
        % test point - soft detection I/Q
        % -----------------------------------------------------------------
        s_pre_det_real = repmat(real(s_pre_det),1,alpha_i_size);
%         s_pre_det_real_soft_pos = exp(-1*(s_pre_det_real-1).^2/(var_dist*xtxPower_mean));
%         s_pre_det_real_soft_neg = exp(-1*(s_pre_det_real+1).^2/(var_dist*xtxPower_mean));
%         s_mean_real    = (s_pre_det_real_soft_pos - s_pre_det_real_soft_neg)./(s_pre_det_real_soft_pos + s_pre_det_real_soft_neg);
        s_pre_det_real_exp = exp(-1*(s_pre_det_real - alpha_i_mat).^2/(MSE(m)*symb_pow));
        s_mean_real        = (s_pre_det_real_exp * alpha_i.')./(sum(s_pre_det_real_exp,2));
        
        s_pre_det_imag = repmat(imag(s_pre_det),1,alpha_q_size);
%         s_pre_det_imag_soft_pos = exp(-1*(s_pre_det_imag-1).^2/(var_dist*xtxPower_mean));
%         s_pre_det_imag_soft_neg = exp(-1*(s_pre_det_imag+1).^2/(var_dist*xtxPower_mean));
%         s_mean_imag    = (s_pre_det_imag_soft_pos - s_pre_det_imag_soft_neg)./(s_pre_det_imag_soft_pos + s_pre_det_imag_soft_neg);
        s_pre_det_imag_exp = exp(-1*(s_pre_det_imag - alpha_q_mat).^2/(MSE(m)*symb_pow));
        s_mean_imag        = (s_pre_det_imag_exp * alpha_q.')./(sum(s_pre_det_imag_exp,2));
        
        s_mean_test    = s_mean_real + 1i*s_mean_imag;
%         real_mean      = mean(s_mean_real);
%         real_var       = mean(s_mean_real.^2 - real_mean.^2);
%         imag_mean      = mean(s_mean_imag);
%         imag_var       = mean(s_mean_imag.^2 - imag_mean.^2);
%         s_mean_var     = real_var + imag_var;
        var_soft1(n) = abs(mean(symb_pow - abs(s_mean_test).^2))/symb_pow;
        %var_soft_test1 = 1 - s_mean_var/xtxPower_mean;
%         var_soft_test1 = 10^(a_cf * var_dist^b_cf + c_cf); % curve fitting version of var estimate
%         s_mean(:,p)    = s_mean_test;
        % -----------------------------------------------------------------
        var_soft(n) = mean(abs(s_mean_test - x_pload).^2)/symb_pow;
        z = demodulate(modem.qamdemod(M_mod),s_mean_test);
        [num]       = symerr(z,x,'column-wise'); % symbol error rate in each frame
        numerr(m) = numerr(m)+num;
        
        waitbar(((m-1)*numTrials+n)/(numTrials*M_LIM))
        
        n = n + 1;
    end
    SER(m) = numerr(m)/(nsymb*numTrials); % SER
    VAR_SOFT(m)  = mean(var_soft);
    VAR_SOFT1(m) = mean(var_soft1);
%     VAR_SOFT2(m) = mean(var_soft2);
    
    m = m + 1;
end
close(h_waitbar)

% End of script.